How can we derive Fourier's Law from quantum mechanics? 
Exact master equation analysis 
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lutionless quantum master equation under three assumptions for the interaction kernel. To second 
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The pure exchange model directly leads to energy diffusion in a weakly coupled spin- 1/2 chain. 
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Introduction. — Fourier's law, connecting the rate of 
heat flow within a body to the temperature profile along 
the flow path, is an empirical law based on observation. 
Despite its fundamental nature, a derivation of this law 
from first principles is still missing [ij. 

In classical systems, extensive numerical simulations 
and rigorous derivations [J, [3, 01 manifested the appli- 
cability of Fourier's law for specifically designed nonlinear 
systems. It is still a great challenge to derive this macro- 
scopic law from microscopic quantum principles. Quan- 
tum simulations of heat flow in spin chains pointed out 
on the validity of Fourier's law, in connection with the 
onset of quantum chaos More recent studies have 
focused on the derivation of this law from Schrodinger 
dynamics [1, S] ■ By using the Hilbert space average 
method, Michel et al. demonstrated the emergence of 
heat diffusion within a quasi degenerate level local Hamil- 
tonian using the truncated Dyson series for the short- 
time displacement operator Q. Using this approach, the 
transition from diffusive to ballistic dynamics at different 
length scales was further explored ^8i]. These works suc- 
cessfully manifested the onset of Fourier's law for a spe- 
cific class of quantum models. The rigorous derivation of 
this law from quantum dynamics for general Hamiltoni- 
ans is still a great challenge. 

In this letter we employ the most general non pertur- 
bative microscopic master equation of the gain-loss form 
to derive quantum-mechanically the Fourier's law. Un- 
like the derivation in Ref. Q, applicable for a class of 
modular designed models, we use the most general local 
Hamiltonian. The diffusive behavior is obtained by mak- 
ing three assumptions on the interacting kernel in the 
microscopic master equation. These assumptions can be 
traced down to the structure of coupling between local 
sites. At weak coupling, going to the long time limit, 
we furthermore demonstrate that the first and second 
assumptions are independent of the specific form of the 
interaction, and the exchange interaction plays the deci- 
sive role in order for the third assumption to hold. 

Master equation. — Consider a system with N lo- 



cal units (particles), each has the same M (could 
be infinite) dimensional eigenspace spanned by bases 
{|n)}, possibly coupled to their own bathes. The 
eigenspace of all N particles is spanned by bases 

{|n) = 1^1)^^2)2 ... |nAf)jv}- The most general 
Hamiltonian of the system and its bath is 



H = Hq + XHj, 



(1) 



where Ho = Hs + Hb and ^fs=SjIi J^m 
is the system Hamiltonian. We assume that the energy 
spectra of the sites are identical, £n(i) = ^nij) > 0, 
and the ground state energy is set as zero. Hb is the 
bath Hamiltonian, consisting degrees of freedom other 
than included in the system. The bath may act locally 
on each site, Hb = Xlfc=i-^B(*)> where HB(i) couples 
to the Hh local site. The second term Hj includes in- 
teractions between system particles, where A character- 
izes the strength of these interactions. It can be gener- 
ally written as Hi = ^ |n) (m| i3„m, where Bnm 's are 
the matrix elements of either system or bath operators. 
The dynamics obeys the Liouville equation ^ ptot (t) = 
-~iX[Hi{t), ptot{t)] = XC{t)ptot{t), where ptot is the den- 
sity matrix of the total system, and we work in the inter- 
action representation Hi(t) = e^^^'iJ/e"*^"*. The Liou- 
ville superoperator C{t) is defined by this equation. The 
super projection operation of interest, p{t) = Vptotit), 
defines the relevant part of the total density matrix for 
the open system. This part exactly satisfies the time- 
local master equation Q 



^p(t)=/C(t)p(i), 



(2) 



given that the initial state is in the relevant subspace, 
VptotiO) = Ptot{0)- The time-convolutionless generator 
IC{t) is in general an extremely complicated object, cal- 
culated using perturbative expansions 0]. Though the 
time-local master equation ([2]) is less well known than 
the Nakajima-Zwanzig equation lull, it is easy to 
show that these forms are equivalent [ij]. In order to 
project the diagonal part of the total density ptot{t) we 
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use the following projection 11, l2| 

Vptot = ^ |n) (n| tr{|n) (n| ptot} pb, (3) 

where pB is the bath thermal equilibrium density matrix, 
and the trace is taken over both system and bath states. 
Note that tr{|n) (n| ptot} = Pn defines the probability 
to find the system in state |n). The time-local master 
equation can be exactly written as [l3| 



dP, 



(4) 



where the rates, or interacting kernels Wnm{t), are com- 
plicated functionals of the interaction Hi{t). Assu ming 
that VC{t)V = 0, which is true for closed systems |11 | 
and for many open systems [9], to second order in Hi 
the expansion lC{t) = /g dsVC{t)C{s)V holds. One then 
obtains 



Vl^mn(i) = 2A2Re j drD^ 
Jo 



>(r)e' 



(Bn-£m)T 



(5) 



nme and L'nm(T) = 



where B„m(T) = e'"^^ B^rne^'"''^ 
trB{Bnin{r)B^nPB}- The total energy En = 
J2iLi eigenvalue of Hs- For a closed system 

Dnm{T) = |i3mn|^, whcrc i?nm 's are c-numbers. 

The gain-loss master equation © was introduced for 
a closed system, and proved vigorously in reference ll|. 
In the long time limit, or Markovian limit, the time- 
dependent interacting kernel W(t) becomes a constant 
matrix for both open or closed systems [l3| ■ In what fol- 
lows we focus on a closed system {Hb is neglected), and 
study energy diffusion between system units due to the 
Hi coupling. The starting point of our derivation is the 
exact equation ([2|) with a non-perturbative kernel W . 

Assumptions and derivation. — First consider a one- 
dimensional system. We employ the nearest neighbor 
interaction form Hj = X^i^i^ ^ih i + 1) where the sym- 
metry V{i,i + 1) ~ V^(i-|-l,i) holds. Owe first assumption 
is that the interacting kernel matrix W takes the same 
symmetry as the interaction, i.e., W = W{i, i + 1) 

and W{i,i + 1) = W{i + l,i). This "localization" as- 
sumption implies that the many-site correlation (kernel 
W) is given by the sum of two-site correlations W{i, i+l). 
This assumption is not trivial as the interacting kernels 
are not lineally related to the interaction Hj. The matrix 
elements of W are therefore given by 

N-l 

{i,i+i;t) n ^ 



i=l 



(6) 

As diagonal elements (n = m) do not contribute to equa- 
tion ^ , they are allowed to be exceptions of the assump- 
tion. Our second assumption describes energy conserva- 



tion between initial and final system states, 

U'n,n,+i(i,i + ni = nii+i and n^+i = rrii , , 

0; Otherwise. 

Besides energy conservation, this condition also implies 
that the local spectra are anharmonic, see discussion af- 
ter Eq. (I17p . In the third assumption, we assume that 
WuiUi+iihi + 1;0 = w{i,t), independent of the rii and 
jT-i+i quantum numbers. The probability to find the *th 
particle in state \ni)^ is 

P„^{i) ^tT{\n,)^{n,\ptot} = J2 

rij ^rii 

Incorporating Eqs. ^ and ^ into ([4]) we obtain 



dt 



= w{i,t)[Pn{i + l)-Pn{i)]- 

w{l-l,t)[P,,{l)-Pr,{l-l)], (9) 



where for convenience we use the short notation P„(j). In 
deriving ([9]) we have also utilized the symmetry W{i,i + 
1) = W{i + 1, i) which holds exactly for a closed system 
at second order [see text after Eq. (O] . Next we write an 
equation of motion for the internal energy at each site, 



du{i) 
dt 



w{i, t){u{i + 1) — u{i)) 
-w{i-l,t){u{i)~u{i-l)). (10) 



The continuous version of this equation is 
du(xA) ■? d f , .du(x,t) 



(11) 



where a is the distance between neighboring sites. Gener- 
alization of this derivation to a three dimensional simple 
cubic lattice is straightforward, leading to Eq. ([TT|) with 
r and V replacing ^ and x. Applying the continuity 



equation for the energy density. 



du{r,t) 
dt 



-V- J(r,i), J 



is the heat current, we exactly obtain the Fourier's law 

J(r, t) = -w{r, t)a'^Vu{r, t) = -K(r, t)VT{v, t). (12) 

Here T(r, t) denotes the temperature profile and k{y, t) = 
w(r, t)a'^Cv with Cy = du/dT as the specific heat. The 
heat conductivity n is essentially time-dependent due to 
the explicit time dependence of the microscopic rates 
w{r,t). In the long time, or Markovian limit, the mi- 
croscopic rates and the temperature profile become con- 
stants, leading to a time-independent relation. Equation 
(fT2|) is the main result of our paper. We emphasize that 
it was derived from the exact master equation Q for a 
generic local Hamiltonian. 
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An exact model in the second order. — We present next 
a model Hamiltonian that exactly satisfies the three as- 
sumption leading to Fourier's law in the Markovian limit. 
We assume a closed system {Hb is neglected), and em- 
ploy a pure exchange interaction form 

^ J,,,+if(*'^+i). (13) 

i=l 

Here is the permutation operator and Ji^i+i are 

nearest-neighbor coupling constants (superexchange for 
spin system), taken as constants Ji.i+i — J and set to 
one in the following discussion. In the two-level case 
(M = 2), since £^'''+^'> = i(l + a, ■ a,+i), Eq. (USD 
is simply the Heisenberg spin-1/2 exchange interaction 
model. The correlation function in ([5]) is given by 

when n ^ m, satisfying the first assumption. Further- 
more, the second assumption is fulfilled here exactly 
as the matrix W has the exact form of Eq. ([7]) with 
'w{i, t) = 2A^t. In order to check the validity of the sec- 
ond order approximation, we compare this result with 
the exact solution for a three unit systems, = 3 and 
M = 2. Let the initial state be [1)^ |0)2 10)3. We can 
get an exact solution for the dynamics for the interac- 
tion (ini), dPi{2)/dt = w{t){Pi{l) + Pi(3) - 2Pi(2)), 
where w(t) — \/2A Mn(2\/2 At) ^ -p^^ weak coupfing (small 

^ ' 2 cos^ (V2At) — 1 ^ ^ 

A), w{t) = 2A^t, which is the same as the second order 
result. While this model satisfies all three assumptions, 
in the long time limit it does not lead to Fourier's law. 

However, the energy difference between spin states 
of site i may slightly differ from that of site i + 1, 

Am^m+iii) = hiiii) - en.+i(i)] - l^niii + 1) - en,+i(« + 
1)] ^ 0. This effect may originate from thermal fluctu- 
ations due to the existence of local bathes at each site, 
which do not contribute to the, diffusive behavior. In 
this case, w„ini+i(«,i) = 2A — ^ " ^-^ — . At short 
times, if the difference A„..„;^j(z) is small, one again 
finds w(t) oc t. In contrast, at long times Wmm+iii) — > 
2TrX^6{Ani,ni^i{i))- When the energy spectrum of each 
spin state forms a band that is dense enough [l^, the 
Fermi's Golden rule is obtained, Wmm+iii) = 27rA^r(0), 
where r(0) is the density of states at zero detuning. 

Thus, this simple exchange model exactly reproduces 
the long time Fourier's law with k = 27rA^a^Cyr(0). 
This result manifests that the local sites do not need 
to acquire exactly identical spectra. It also provides a 
microscopic explanation of the validity of Fourier's law 
in spin chains with Heisenberg-type interactions [l^ . 

Long Time (Markovian) Limit. — Next we show that 
at second order in the interacting kernel the first two 
assumptions [Eqs. ©-([T])] are model-independent in the 



long time limit, and the third assumption is valid for fam- 
ily of interactions with exchange effect. We consider the 
most general interaction V{i, i + 1) for a one-dimensional 
system. It is easy to show that 

In the long time (Markovian) limit Eq. ([5]) reduces to 

W^m„(t ^ ^) = 27rA2 5{E^ - E^), (16) 

using '^'"[^"Tgf")^*"' T^^iEn — -Em) when i — > oo. As- 
suming the bound states are non-degenerate, a somewhat 
tedious calculation shows that the first assumption is sat- 
isfied, and the off-diagonal matrix elements are given by 

^^TiiUi-^-i ^mirriij^i (^i ^ ^" ^) 

(z,z-f l)|'5(A„^„^^^,„^„^^J, 

(17) 

where A„^„.^^ = e„^ + 6^+^ - - Cmi+i- Since 
e„ > 0, the delta function in Eq. (fT7|) implies that 
^rii = Cmi+i and erii^-i = ejTii, provided that the en- 
ergy spectra are anharmonic. Therefore, in the non- 
degenerate case, only transitions between the quantum 
states Ui = rui+i and n^+i = rui exist, in accordance 
with the second assumption. In contrast, when the 
energy spectra are strictly harmonic^ the delta func- 
tion can be satisfied for large number of combinations, 
en, - Cm; = em.+i " ^n.+i = j^^o- Hcrc t^o is the en- 
ergy difference between subunit states (equal for all N) 
and j is an integer. In this case the second assumption 
does not hold, and we cannot derive Eq. Q and the 
subsequent result (|12p . Thus, interestingly, in order to 
derive Fourier's law, the system spectrum should be an- 
harmonic, in accordance with classical results 0, 

At second order, the first two assumptions are there- 
fore the results of the Markovian limit. As discussed 
above, we can expect that there is a dense band structure 
around each level e„, originating from bath fluctuations 
or quantum tunneling between sites. Defining the den- 
sity of states at the energy difference around A as r(A), 
we obtain [l3| 

W„.„.+ i(i) = 27rA^r(0) |Kini + i,ni+in, («, * + 1)|^ , (18) 

which manifests that the third assumption is valid, de- 
pending on the form of specific Hamiltonians. We calcu- 
late next the microscopic rates for different types of inter- 
actions. Assuming translational symmetry in the system, 
we need only discuss a pair of sites, for instance between 
site 1 and 2, Kms.nsni (1,2) = (nin2| V{1,2) \n2n1). 

In general, particles at each site may be electrons or 
atoms. To simplify, we consider the case with one particle 
at each site. Our first example is the short-range delta 
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interaction, V(l,2) = gS{ri2), where ri2 = |a;2— 
and Xi{x2) is the coordinate of the first (second) particle. 
Physically, this interaction describes particles that move 
almost independently in the site interior, while collisions, 
leading to energy exchange between particles, occur at 
the edge points xi ^ xi- This picture is appropriate for 
describing phonon collisions in solids 17[. For this type 
of interaction the matrix elements become 



We also calculate the interaction 1^(1,2) for the har- 



monic potential Virx^'i) 



We find that when ajd 



Ki«2,n2«i (l,2)=.g / dx|0„i(x)|^ |0„2(2;)|^ . 



(19) 



This integral is almost independent of quantum numbers 
n\ and n2 for many systems 18|. Therefore, the third as- 
sumption is generally valid for short-range interactions. 
For example, in one-dimensional infinite square well 

(width d) with the wave function (j>n{x) — y^sin^i^, 

the matrix element Vn-^n2in2ni (li2) = l/d is completely 
independent of quantum numbers, assuming d > a so 
that two particles could collide. 

Our second example is the long-range interaction. We 
choose the general form ^(1, 2) = y(ri2)(Ai + (^'2^), 
where the A's are constants and £^^'^^ is the spacious 
exchange operator flQl. Considering the first order con- 
tribution to this potential, 1^(1, 2) = |x2 — a:i|(Ai -I- 
A2^^^'^^), the matrix elements Kiin27n2ni (li2) become 
the sum of the direct {D) and exchange (E) terms, 

^1712 '»2"1 (Ij 2) = 



Ai / dxidx24>l^{xi)(l)l^{x2)\x2 - Xi\(j)nAx2)4>n2{xi), 



nin2 Jn2"l (^i 2) — 



A2 / dxidX2\(j)nA^l)f \'t>7i2{x2)\'^ \X2 - Xi\ . (20) 



In what follows we take Ai = A2 = 1. If the overlap 
between the wave functions of two particles is zero, the 
direct integral diminishes because of the orthogonality 
of two states, and the exchange integral {x2 > xi) is 
Kun2'"2ni (1,2) = a + {n2\x2 1^2) - (nil^i Im), where 
x'2 = X2 — a is the relative coordinate of particle 2. If the 
wavefunction 0„(x) has a well defined parity, {n\ x\n) = 
for the relative coordinate of each particle, implying that 



14,in2jn2ni (1,2) = o is a constant. This satisfies the third 
assumption. Note that the result is independent of the 
details of the local Hamiltonian. 

If the overlap is not zero, yet small, the direct inte- 
gral is expected to be small. Considering again the one- 
dimensional infinite square well, when a/d < 1, two par- 
ticles will overlap. Numerical calculations show that for 
quantum numbers ni , n2 ranging from 1 to 20 the direct 
integrals are less than one tenth of the exchange integrals 
and Kiin2,n2ni (1,2) is almost a constant, very close to a 
(derivation around 3% when 0.4 < a/d < 1). The bigger 
the values of a/d is, the smaller the ratio between the 
direct and the exchange integrals is. 



varies from 0.7 to 0.95 the exchange integral is almost a 
constant with deviations from 14% to 4%, and the ratio 
between the direct and exchange integrals is 0.14-0.1. 

The temperature dependence of the heat conductivity 
K results from the interplay between the specific heat and 
the interaction kernel. The specific heat, defined per unit, 
can be easily calculated for different models, e.g., for a 

spin chain of ujq spacing, Cy = ^e"o/^(e'^«/^ + 1)"^ 
In contrast, our approach does not directly bring in the 
temperature dependence of the kernel IVmn- Physically, 
since Vn-^n2 ,n2ni (1, 2) is a function of the intersite separa- 
tion, the lattice vibrations can modify it, thus introduce 
temperature dependent transition rates. This effect can 
be successfully included by phenomenologically introduc- 
ing a temperature dependent interaction [2^, [2l[ . Here 
we adopt the simple form Vt(''12) = V{ri2)ji'ri2,T), 
where T may be the average temperature of site 1 and 
2, and 7(7-12, = exp(— a/(T)ri2). In the case of in- 
teracting ions [20|, a is inversely proportional to the 
Fermi velocity and /(T) = T. When /(T) 00, 
VT{ri2) af'(T) ^(0)^('"i2)- Therefore, for high temper- 
atures K oc l//(r)^, as the lattice specific heat typically 
saturates. For the ionic solid of Ref. [13] one obtains 
K oc l/T^, in agreement with standard expectations [l7| . 

Summary. — We have presented here a microscopic 
quantum derivation of Fourier's law of heat conduction 
that is not limited to specifically designed models. The 
derivation relays on three assumptions that are satisfied 
in weak coupling and at long times for a family of ex- 
change interaction potentials. Our analysis naturally im- 
plies that energy diffusion cannot emerge in harmonic 
models, in agreement with the behavior of classical sys- 
tems. 
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